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Abstract 

This note discusses a theoretical issue regarding the application of the "Modular Response Analysis" 
method to quasi-steady state (rather than steady-state) data. 



1 Introduction 

The reverse engineering problem in systems biology is, loosely speaking, that of unraveling the web of 
interactions among the components of protein and gene regulatory networks. A major goal is to map out 
the direct functional interactions among components, a problem that is difficult to approach by means 
of standard statistical and machine learning approaches such as clustering into co-expression patterns. 
Information on direct functional interactions throws light upon the possible mechanisms and architecture 
underlying the observed behavior of complex molecular networks. 

An intrinsic difficulty in capturing direct interactions in intact cells by traditional genetic experiments, 
RNA interference, hormones, growth factors, or pharmacological interventions, is that any perturbation 
to a particular gene or signaling component may rapidly propagate throughout the network, thus causing 
global changes which cannot be easily distinguished from direct (local) effects. Thus, one central goal in 
reverse engineering is to use the observed global responses (such as steady-state changes in concentrations 
of active proteins, mRNA levels, or transcription rates) in order to infer the local interactions between 
individual nodes. One potentially very powerful approach to solve the global to local problem is the 
Modular Response Analysis (MRA) or "unraveling" method proposed in [5] and further elaborated upon 
in [g[Tl |3 E] (see [TpJIU for reviews). 

The MRA experimental design compares the steady states that result after performing independent 
perturbations to each "modular component" of a network. These perturbations might be genetic or 
biochemical, or (in eukaryotes) they might be achieved through the down-regulation of protein levels by 
means of RNAi. This latter experimental approach to MRA was the one taken in [7J, which quantified 
positive and negative feedback effects in the Raf/Mek/Erk MAPK network in rat adrenal phcochromo- 
cytoma (PC-12) cells. Using the formulas given in [S] and pQ, the authors of [7J uncovered connectivity 
differences depending on whether the cells are stimulated with epidermal growth factor (EGF) or instead 
with neuronal growth factor (NGF). There are a couple of subtle theoretical gaps, however, when ap- 
plying MRA algorithms to data like that employed in [7] . The main gap is due to the fact that the data 
fed into the MRA algorithm included non-steady state measurements. Specifically, for EGF stimulation, 
network responses were measured at the peak of Erk activity (at 5 minutes) and not at steady state. This 
note fills that gap, providing a theoretical justification for the use of quasi-steady state information. 

1.1 Mathematical formulation 

We assume that there are n quantities Xi (t) that can be in principle measured, such as the levels of activity 
of selected proteins, or the transcription rates of certain genes. These quantities are thought of as state 
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variables in a dynamical system and are collected into a time-dependent vector x(t) — (xi(t), . . . , x n (t)). 
The dynamical system is described by a system of differential equations: 

&i = fi(xi,...,x n) p li ...,p m ) 

%2 = f2(Xl,...,X n ,pi,...,p m ) 
%n In [%1 ) • • • ; %n : Pi j • • • j Pm) 

or, in more convenient vector form, 

x = f(x,p) 

(dot indicates time derivative, and arguments t are omitted when clear). The p^s are parameters, 
collected into a vector p = (pi, . . . ,p m ). These parameters can be manipulated, but, once changed, they 
remain constant for the duration of the experiment. An example would be that in which the variables Xi 
correspond to the levels of protein products coresponding to n genes in a network, and the parameters 
reflect translation rates, controlled by RNAi. Another example would be the total levels of proteins 
whose half-lives are long compared to the time scale of the processes being described by the differential 
equations, such as phosphorylation modifications of these proteins in a signaling pathway. 

The ultimate goal is to obtain, for each pair of variables Xi and Xj, the relative signs and magnitudes 
of the partial derivatives 

dxj ' 

which quantify the direct effects of each variable Xj on each variable Xi. 

The critical assumption for MRA, and indeed the main point of [5] 6, 9J, is that, while one may not 
know the detailed form of the vector field /, often one does know which parameters Pj directly affect 
which variables x t . For example, Xi may be the level of activity of a particular protein, and pi might 
be the total amount (active plus inactive) of that particular protein; in that case, we know that pi only 
directly affects X{. 

In the standard version of MRA, one first measures a steady state x corresponding to a "wild type" 
vector of parameters p, that is f(x,p) = 0. Subsequent perturbations are separately performed to each 
entry of p 7 and a new steady state is measured, one for each such perturbation. Using these data (and 
assuming a that certain independence condition which we review later is satisfied), it is possible to 
calculate, at least in the ideal noise-free case, the Jacobian of /, evaluated at (x,p), up to a scalar 
multiplicative factor uncertainty on each row. (Such uncertainty is unavoidable when using only steady 
state measurements, since multiplying a row of the vector field / by a nonzero constant does not affect 
the location of steady states.) A variation of MRA is possible, which allows for the use of non-steady 
state, time-series data. However, this alternative method, developed in [9J, requires one to compute time 
derivatives, and hence is hard to apply when time measurements are spaced far apart and/or are noisy. 
An intermediate possibility is to use quasi-steady state data, meaning, just as in the experimental setup 
of [7], that one employs data collected at times when a variable has been observed to attain a local 
maximum or local minimum. That is the case addressed in this note. 

More precisely, we will consider the following scenario. For any fixed variable, let us say the ith 
component Xi ofx, we consider some time instant U at which Xi(t) is zero. Under the same independence 
hypothesis as in the classical MRA case, plus the nondegeneracy assumption that the second time 
derivative Xi(ii) is not zero (so that we have a true local minimum or local maximum, but not an 
inflection point), we show here that the MRA approach applies in exactly the same manner as in the 
steady-state case. Specifically, the ith row of the Jacobian of /, evaluated at the vector (x,p), is recovered 
up to a constant multiple, where x = x(ti) is the full state x at time ii. The main difference with the 
steady-state case is that different rows of / are estimated at different pairs (x,p), since the considered 
times U at which each individual ii(t) vanishes are in general different for different indices i, and so the 
state x is different for different i's. 
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2 Using quasi-steady state data 



We fix an index i 6 {1, . . . , n}, and an initial condition x(0), and assume that the solution x(t) with this 
initial condition and a given parameter vector p has the property that, for some time i = i,-, we have 
that both ii(t) — and Xi(i) ^ 0. At the instant t = i, Xi achieves a local minimum or a local maximum 
as a function of t. We describe the reconstruction of the ith row of the Jacobian of /, which is the same 
as the gradient V/,, where fi is the ith coordinate of /, evaluated at x — x and p — p, where x = x(t). 

To emphasize the dependence of the solution on the parameters (the initial condition x(0) will remain 
fixed), we will denote the solution of the differential equation x = f(x,p) by x{t,p). The function x(t,p) 
is jointly continuously diffcrcntiablc in x and p, if the vector field / is continuously differentiable (see 
e.g. [8], Appendix C). Note that, with this notation, the left-hand side of the differential equation can 
also be written as dx/dt, and that x(i,p) = x. 

We introduce the following function: 

dx 

Note that a(i,p) = 0. Also, 

doc O^X' 

— (t, P ) = -Qjrfop) = VfMt,p),p)f(x{t,p),p). 

The assumption that Xi(i) ^ when p = p means that ^(?, P) 0. Therefore, we may apply the 
implicit function theorem and conclude the existence of a mapping r, defined on a neighborhood of p, 
with the property that 

a(r{p),p) = for all p « p 
and t(p) = t (and, in fact, t = r(p) is the unique value of t near t such that (dxi/dt)(t,p) — a(t,p) = 0). 
We define, also in a neighborhood of p, the differentiable function 

<f(p) = x(r(p),p) 

and note that ip(p) = x. Observe that, from the definition of a, we have: 

fi{v{p),P) = for all p w p. (1) 



We next discuss how to reconstruct V/,(S,p), up to a constant multiple, under the assumption (as 
in [S]) that it is possible to apply n—1 independent parameter perturbations to all species different from 
the ith one. This discussion is basically identical to that for the steady state case, given in [5j[l][2]. 

Mathematically, we assume that there are n — 1 indices ■ ■ ■ ,jn-i with the properties that (a) 

fi does not depend directly on any pf dfi/dpj = 0, for j £ {ji,j2, • ■ • ,jn-i}, and (b) the vectors 
Vj = (dip/dpj)(p), for these j's, are linearly independent. Assumption (a) is structural, and is key to 
the method and nontrivial, but assumption (b) is a weak genericity assumption. 

We then have, taking total derivatives in ([1]): 

Vf l (x,p)v j = 0, j e {p n ,p h , ■ ■ ■ ,Pj„^A- 

Thus, the vector V/i(x,p) is orthogonal to the n—1 dimensional subspace spanned by {v\, . . . , 
and hence is uniquely determined up to multiplication by a positive scalar. Another way to phrase this 
is to say that \7fi(x,p) is in the (one-dimensional) left nullspace of the matrix A whose rows are the 
i>i's, or that (if nonzero) the transpose of this gradient can be found as an (any) eigenvector associated 
to the zero eigenvalue of the transpose of A. 
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Numerical approximation by finite differences 



Approximating the vectors Vj by finite differences, one has that Vfi(x,p) is approximately orthogonal 
to these differences as well. Explicitly, suppose that we approximate Vj = (d(<p/dpj)(p) by: 



where h is small and where ej is the vector having a one in the jth position and zeros elsewhere. Then, 
V/j(5,p) is (approximately) orthogonal to the differences 



which form a set of n — 1 linearly independent vectors (if h is small). A simple matrix inversion (after 
fixing an arbitrary value for one of its entries) allows the computation of V/i (x, p) . Observe that division 
by the potentially small number h is not required in performing these operations, In fact, no knowledge 
whatsoever about parameter values is needed by the algorithm. 

Note that ip(p + hej) is the state x(t) at the time t at which the particular coordinate Xi achieves a 
local extremum value, if the parameters have been perturbed to p — p + hej. To be more precise, t is 
the unique time close to t such that ii(t) = when parameter vector p is being used. Theoretically, we 
must have p w p, so h must be very small, but, in practice, quite large perturbations of p also work fine. 

3 A simple numerical example 

We illustrate the calculations with a very simple example, the following system (writing x instead of Xi 
and y instead of X2)' 



with initial state (0,0) and reference parameter p = 2. This might represent the simplified dynamics of 
a two-gene network, in which the first gene enhances the expression of the second gene, which in turn 
represses the rate of expression of the first one, there is a constitutive rate of production of the second 
gene, and both protein products decay at rate 3 sec -1 . The single parameter p may represent a promoter 
strength, and we assume that there is a way to perturb it (perhaps by duplication or sequence change) . 
The solid lines in Figure Q] (and also in Figures and [3]) show plots of the solution coordinates x{t) and 



Let us pose the following problem: not knowing the above equations, estimate the relative strength 
of the second gene's effect on the rate of expression of the first one. The only data to be used are the 
levels of both gene products (x(t) and y(t)) at the time when x(t) achieves its local maximum. We do 
assume known the fact that the parameter p affects directly only the rate of expression of the second 
gene, not the first. Observe that the maximum of x is attained at i k 0.5275, and the values there 
are (approximately) x(t) = 1.6553 and y(t) = 1.0f38. The gradient V/i of — 3x + j^, evaluated at 
(1.407, 1.3695), has the true (but unknown to the algorithm) value: 



Next, we perform the "experiment" in which p is up-perturbed by 25%. With the new parameter 
p = 2.5, we obtain plots as shown by the dotted lines in Figure [TJ Now the maximum of x is attained 
at t w 0.4268, and the values there are x{t) = 1 .407 and y{t) = 1 .3695. Letting 8 = (f .407, f .3695) - 
(1.6553, f .Of 38), the unknown (to the algorithm) gradient V/i is known to be (approximately) orthogonal 



- [ip{p + hej) - x) , 



ip{p + hej) - x, 



x = — 3x H 

l + y 

y = px + 1 — 3y 



y(t)- 
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Figure 1: Trajectories; dashed is perturbed motion with 25% change in parameter 

to S. Any vector perpendicular to S must be a multiple of (—3, —2.3455). (We normalized the first entry 
to -3 merely in order to compare our result to the true gradient; the algorithm does not know the value 
"-3". In practice, however, one may assume that the first entry of the vector is negative, reflecting 
degradation or dilution effects, so the algorithm will give the correct sign for the second term, as well as 
its magnitude relative to the rate of degradation or dilution.) The relative error in our estimate is less 
than 5%. 

Even larger perturbations may be performed. For example, a 50% perturbation from p — 2 to p = 3, 
provides the dashed lines in Figure [D Now the maximum for x is attained at t w 0.4658, and there 
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Figure 2: Trajectories; dashed is perturbed motion with 50% change in parameter 

x(t) = 1. 5103 and y(t) — 1.2073. The estimated gradient is now (—3, —2.2476), which gives a relative 
error of less than 9%. Finally, a 100% perturbation to p = 4 provides the dashed lines in Figure [3] Now 
the maximum for x is attained at t « 0.4268, and there x(t) = 1.4071 and y(t) = 1.3695. The estimated 
gradient is now (—3, —2.0936), which gives a relative error of about 15%. 

Remarks 

As its name implies, one of the main advantages of the MRA method in the steady-state case is that 
only "communicating intermediates" in-between "modules" need to be measured (for example, just the 
active forms of Erkl/2, Mekl/2 and Raf-1, in [7]). Here, we only carried out the analysis in the case in 
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Figure 3: Trajectories; dashed is perturbed motion with 100% change in parameter 



which all the variables x\ can be measured. In the general case, if one assumes that hidden (internal) 
variables are at quasi-steady state at the same times as the communicating variables, then an implicit 
function argument as in [6] allows one to reduce to the present situation, by writing the hidden variables 
in terms of the communicating quantities. However, there is no reason for the method to work when the 
hidden variables do not have this property. 

Also, we assume perfect "noise free" data. The analysis of noise performed in [1] carries over with 
no changes to the quasi-steady state case. 
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